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Abstract 

Kinetic energy functionals of the electronic density are used to model large systems in the context 
of density functional theory, without the need to obtain electronic wavefunctions. We discuss the 
problems associated with the application of widely used kinetic energy functionals to non-periodic 
systems. We develop a method that circumvents this difficulty and allows the kinetic energy to be 
evaluated entirely in real space. We demonstrate that the method is efficient [O(N)] and accurate 
by comparing the results of our real-space formulation to calculations performed in reciprocal space, 
and to calculations using traditional approaches based on electronic states. 
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Electronic structure calculations based on density functional theory]!], (DFT) and em- 
ploying an approximate kinetic energy functional || f|, [5], |6|, [7], ||, |9] have been shown to 
yield accurate energetics for a variety of physical systems, with a considerably smaller ex- 
penditure of computer effort than traditional schemes. A drawback of many existing kinetic 
energy functionals is the necessity to evaluate them in part in reciprocal space. The benefits 
of doing electronic structure calculations exclusively in real space are several. Foremost 
among them is their ability to simulate efficiently isolated systems, whereas reciprocal space 
methods would require a large supercell with a significant portion of the volume devoted to 
uninteresting vacuum. A direct extension of this feature is the possibility of using arbitrary 
boundary conditions rather than the strict periodic boundary conditions underlying recip- 
rocal space approaches; this should be of paramount importance in dealing with complex 
structures which cannot be accommodated by simple periodic boundary conditions, such as 
dislocations, cracks, etc. Finally, real space methods can be readily parallelized for efficient 
computations on parallel computer architectures. 

In this article we examine the reasons for the evaluation in reciprocal space of many kinetic 
energy functionals and propose a new functional form that can reproduce the energetics of 
those functionals but does not require any reciprocal space evaluations. The performance of 
the new functionals is evaluated and compared to the existing reciprocal space functionals 
as well as to the traditional approaches based on the calculation of electronic states. 

DFT in its usual guise consists of solving the Kohn-Sham (KS) equations^. The fictitious 
KS wavefunctions allow the exact evaluation of T s [p], the kinetic energy of non- interacting 
fermions with a density p(r). Recent calculations have shown that for certain systems the 
solution of the KS equations, usually the most computationally demanding part of electronic 
structure calculations, can be bypassed using an approximate form for T s [p] ||. A fair number 
of approximate non-interacting kinetic energy functionals have been proposed for use within 
such orbital- free methods. A subset of them are similar in form|4], |5j, |], [7], ||, sharing the 
following common traits: 

1. A major ingredient of T s [p] is the Thomas-Fermi (TF) energy JIT]]: 



where Ctf — jq (3vt 2 ) . Here and throughout this article atomic units (h = m e = 
e = 1) are employed. The TF energy approximates the kinetic energy in an element of 
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space dr with that of a homogeneous non- interacting electron gas with a density p(r). 
Hence for homogeneous densities p(r) = p , Ttf\p(t)] is exact. 



2. Another important contribution is the von Weizsacker energy given by: 

T vW [p] = ~\f vW) V 2 yW) dr (2) 

It can be readily shown that the von Weizsacker energy yields the correct kinetic 
energy for a non-interacting fermion density that consists of a single orbital, i.e. a 
one- or two-electron density. Also, for any density, the von Weizsacker energy yields 
the energy that a system of non-interacting bosons of density p(r) would have. This 
term is a lower bound of T s [p] (K 



3. The response of a homogeneous non-interacting Fermi gas to a small perturbation is 
known exact lyfll2"|. This response can be related to the second functional derivative of 
T s [p] evaluated at uniform density p(r) = p . The Fourier transform of this functional 
derivative is given by: 
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where XLind(q) is the Lindhard response function, jF[/(r)] = f f(r)e lk ' r dr denotes the 
action of the Fourier transform on the function /(r), q = k/2/jp, and q — |q|. In order 
to satisfy Eq.(||) in addition to items [I] and above, the kinetic energy functionals 
under consideration here include a term T#-[p], which will be hereafter referred to as 
the kernel energy: 



Tk[ P ] = / f(p(r))K(\r-r'\)g(p(r'))drdr' 



(5) 



where for the moment f(p), g(p) are arbitrary functions that can be chosen to satisfy 
known limits of the exact T s [p\. 

The total kinetic energy functional is taken to be the sum of these terms: 



T s [p}~T TF [p}+T vW [p}+T K [p} 



(6) 
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By plugging Eq.@ into Eq.(|3|), it is seen that T s [p] exhibits the correct linear response, 
provided: 
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The functional derivatives of T FF and T v w are well known, and the second functional deriva- 
tive of Tk can, by design, be easily evaluated, so that Eq. (^) takes the form: 
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For any choice of the functions / and g, K(q) can be readily chosen such that the total 
kinetic energy functional exhibits the correct linear response. The different kinetic energy 
functionals considered in this paper differ mostly in the choice of the functions / and g. 

In order to use kinetic energy functionals in actual electronic structure calculations, an 
algorithm must be developed for their action on discrete representations of the charge density. 
It is clear that Tt F [ P (t)] and T„y^[p(r)] can be computed easily and efficiently on a grid in 
real space: 



T TF [{ Pl }] = nc TF Y,p' 
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where Q is the volume per grid point, { Pi } denotes a discrete representation of p(r) on a grid 
in real space, and Ay is a discrete representation of the Laplacian operator. In principle, 
we could also compute the kernel energy T^[p(r)] in real space as follows: 



t k [{ Pi }] = n 2 J2f (Pi) Kite - TjDg ( Pj ) 



(12) 



However, this is impractical because of the specific form of K(r). K(r) is the Fourier trans- 
form of K(q) (see Eq.(Bj)), which is a non-trivial Fourier transform. Herringjr3| has shown 
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how it may be evaluated numerically. The quantity K(r) is shown in Fig. p], indicating 
that K(r) does not decay rapidly. In order to evaluate 



accurately, the integral needs to be evaluated over a sphere centered at r of a radius equal to 
A, with A a large value. Such a convolution integral can be computed efficiently in reciprocal 
space: 



Starting from K (r) and h(r) in real space, three Fast Fourier Transforms (FFTs) are required 
to evaluate the convolution: two forward transformations, and one reverse transformation. 
The drawback of computing the kinetic energy functional with FFTs is that this approach 
maps the problem to a periodic tiling of the system of interest. This periodicity can have 
consequences on the resulting physics that in many cases are undesirable and can lead to 
erroneous results. 

The problem with evaluating the kinetic energy functionals described above in real space 
lies in the efficient evaluation of convolution integrals. Convolution integrals with a long- 
ranged kernel can sometimes be evaluated efficiently in real space in an indirect way. For 
instance, the electrostatic potential of a charge distribution p(r) is the convolution of p(r) 
with the very long-ranged 1/r: 



The integral in this expression can be computed efficiently in real space by solving the 
Poisson equation: 
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In general, as an alternative to evaluating the convolution 
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one may equivalently solve the integral equation: 




(18) 
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In discretized form, this integral equation is a standard linear problem: 

HiM = fi (19) 



which can be solved by an iterative linear solver |TJj]; the efficiency of such approaches depends 
on the efficiency of multiplying an arbitrary vector by a matrix. For the case of the Poisson 
equation, the matrix Hij becomes a discrete representation of the well-localized Laplacian 
operator, and hence matrix- vector multiplications are efficient. However, when the integral 
equation that corresponds to the convolution of Eq. fl5|) is constructed, the resulting function 
Hk{t) is found to be long-ranged. The long-ranged nature of both K(r) and Hk{t) is 
due to the logarithmic divergence of the slope of Xund{o) at g = 1, which causes long- 
ranged oscillations to appear in its Fourier transform. Thus, matrix-vector multiplications 
by = Hx(\rj — Tj\), and the solution of the integral Eq.(|T8D, are inefficient. 

We are therefore interested in developing a method that circumvents these difficulties. 
As a first step toward this goal, we note that the kernel appearing in the class of kinetic 
energy functionals under consideration, K(q), can be fit well by the rational function: 

Zt i N ^ 2 ± ■ ■ ■ ± N ^ 2m torn 
= A, + A«» + - + JW <20) 
with appropriate choices of the real coefficients iVj and Di. The odd powers of q are omitted 

because in the Taylor expansions of K(q) about and oo only even powers of q appear. The 

quality of the fit is shown in Fig. ^ for several values of m. Next we note that K(q) can be 

separated into terms of the form: 

m p 2 

(21) 

where the Pj and Qj are now complex numbers. With this expression, the convolution of a 
function /(r) with K(r) in reciprocal space becomes: 

V(q) = K{q)f(q) = V 1 (q) + --- + V m (q), (22) 

where 



q 2 + Qj 

The V%(t) can be computed efficiently in real space by solving 
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that is, Vj(r) is the solution of a complex Helmholtz equation. A shortcut in computing 
the Vj(r) results from the fact that K(q), and thus the sum of the Vj(r), is purely real. 
For every pair of coefficients {Pj,Qj}, another pair {Pk,Qk} = {Pj,Q*j} must also appear 
in the expansion. It follows that Vfe(r) = V*(r), and thus only half of the Vj(r) need be 
computed. A generalization of the kinetic energy functional with the form of Eq.(^) has 
been developed by Wang et al.[[J. These functional can also be treated in real space with 
the present method, as discussed in the Appendix. 

The following issue regarding the form of the approximation of Eq. ( P0|) deserves further 
discussion: an important feature of the Lindhard response function lies in the logarith- 
mic singularity in its slope at q = 1. As discussed above, this singularity manifests itself 
mathematically in the long-ranged nature of K(r) and Hx(r). Not surprisingly, this singu- 
larity also has important physical consequences, such as Friedel oscillations and the Kohn 
effect [p~5f| . The approximate kernel K(q) does not exhibit the singularity. It may seem then 
that from a physical standpoint, K(q) may not adequately describe the kinetic energy of the 
electron gas. However, it should be noted that in a discrete representation of the problem, 
the exact singularity at q = 1 will not be seen. Furthermore, at non-zero electronic temper- 
atures, however small, the singularity in Xiind{q) disappears. Thus, one could think of the 
fitting form K(q) as representing K(q) for a small but finite electronic temperature. The use 
of a fictitious, finite electronic temperature is a trick routinely employed to aid numerical 
convergence in standard DFT calculations of metals. 

We have therefore reduced the problem to solving the complex Helmholtz equation, which 
by itself provides a special challenge due to the fact that the operator is non- Hermit ian. 
Typical iterative methods for solving linear systems, like the conjugate- gradient algorithm, 
fail for non-Hermitian matrices. The complex Helmholtz equation is an important problem, 
arising frequently in the context of electrodynamics. Several iterative methods have been 
developed for the special class of complex symmetric matrices, into which complex Helmholtz 



operators fall|T6|, [L7|]. For the present tests, the biconjugate- gradient algorithm, specialized 
to complex symmetric matrices |16[], has been employed to solve Eq. (p4|) . Each iteration of 
this method requires an amount of computation that scales linearly with the system size; 
thus, if the number of iterations required to converge a solution of the Helmholtz equation 
does not vary significantly with the grid size, the entire method for calculating the total 
energy scales linearly [O(N)] with the system size (N). 
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Solving the discretized version of Eq. ([24]) in real space suffers from another source of 
inaccuracy, beyond that introduced by the fitting of K(q) by K(q). A discretized version 
of the Laplacian operator must be employed, which, in reciprocal space, deviates from the 
exact q 2 behavior. The two sources of error can be easily separated for the purposes of 
numerical tests. To measure the error due to the approximation of K(q) by K(q) alone, the 
kernel energy Tjc [p] can be computed with the reciprocal space convolution method, Eq. (|I4]) , 
but using K(q) instead of K(q). The error due to the use of the discretized Laplacian, in 
addition to the fitting error, is present when the full real space evaluation method is used. 
We found that the error due to a fourth-order discrete Laplacian operator was negligible 
compared to the error introduced by the fitting of K, and throughout the numerical tests 
this fourth-order Laplacian was employed. The kinetic energy evaluated with the present 
real-space method is denoted by T^[p], while the kinetic energy evaluated in reciprocal 
space with K(q) is denoted by T K [p\. The kinetic energies T K [p] and Tj>[p] are evaluated 
and compared for a realistic set of charge densities p(r) in Fig. |3|. The Tk\p\ used is due to 
Wang et al.0, and has parameters {a, f3} — | ± and for T^[p], successive fitting orders 
m = 2, 3, and 4 were tried. The fitting coefficients for m = 4 are given in Table |. The 
charge densities considered are generated by minimizing the total energy: 

Etot = T TF [p] + T vW [p] + T K [p] + E H [p\ + 

E lon [p]+E xc [p] (25) 

where Eh-, Exc, and Ei on are the Hartree, exchange-correlation, and electron-ion interaction 
terms, for a bulk fee aluminum system, with a wide range of lattice constants. Aluminum 
was represented by the Goodwin-Needs- Heine local pseudopotential JT%|] , and exchange and 



correlation were treated with the LDA[|T^]. At each lattice constant, after minimizing the 
electronic energy with the kernel energy represented by Tk[p\, the kernel energy is also 
computed with T^[p] and compared to Tx[p]- As can be seen from Fig. H, the use of more 
terms in the fitting (higher m) results in a smaller deviation, and at m = 4 the accuracy is 
quite satisfactory (less than 1% deviation for the entire range of lattice constants considered.) 

We have performed a different set of tests, in which at each value of the lattice constant 
the total energy with the kernel energy represented by either T^[p] or T^[p] is minimized 
with respect to the density. In the test discussed earlier the density was fixed to that 
obtained from minimization of the total energy using Tk[p\- The present test also differs 
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from the last in that the error of the discretized Laplacian is also present in the evaluation 
of TWfp] in the real space calculations. Thus in this calculation, total energy calculations 
are performed fully and self- consistently in real space, and are compared with reciprocal 
space results. The equilibrium lattice constant and the bulk modulus obtained from these 
calculations are given in Table [IT]. As seen from this Table, m = 4 offers an acceptable level 
of accuracy: the values obtained with this approximation to the true kernel are essentially 
the same as those from the exact kernel, differing only by 0.1% for the lattice constant and 
by 0.5% for the bulk modulus. 

In conclusion, the convolution integrals that appear in the class of kinetic energy function- 
al under consideration in this paper cannot be efficiently evaluated directly in real space, 
and the corresponding inverse integral equations cannot be efficiently solved. The kernels 
of the convolutions can be approximated by a sum of sub-kernels. This approach makes 
it possible to evaluate the convolution efficiently in real space by solving the inverse inte- 
gral equations that correspond to the sub-kernels, which are complex Helmholtz equations. 
We have demonstrated that this method yields excellent results in numerical tests, in the 
sense that the error introduced by the real space method is negligible compared to the error 
inherent in the approximate kinetic energy functionals. 
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APPENDIX A: REAL-SPACE EVALUATION OF KINETIC ENERGY FUNC- 
TIONALS WITH DENSITY-DEPENDENT KERNELS 

Wang, Govind, and Carter |J (WGC in the following) have developed a class of kinetic 
energy functionals that include the Thomas-Fermi and von Weizsacker terms, as well as a 
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term analogous to Eq.(|5]), but with a density-dependent kernel in the convolution: 



T k [p]=CtfX 
f (p(r)) K(p(r), p(r'); |r - r'\)g(p(r'))drdr' 



(Al) 



where 



K(p(r), p(r'); |r - r'|) = K(p(r% p(r); |r - r'|). 



Numerical tests indicate that these kinetic energy functionals are more transferable to sys- 
tems that deviate significantly from the bulk, like surfaces. 

Even by utilizing FFTs, a straightforward evaluation of the convolution in Eq. (|A1| ) 
would require 0(N 2 ) operations, where N is proportional to the size of the system. WGC 
have demonstrated how this convolution can be efficiently, but approximately, evaluated. 
By Taylor expanding K(p(r), p(r'); |r — r'|) with respect to p(r) about some chosen average 
density p, one obtains: 

iir(p(r),p(r');|r-r'|) = J K'o(|r-r , |) 
+ K 1 (\r-T f \)[Ap(T) + Ap(i>)] 
+ ^ndr-r'DfA^W + Ap 2 ^ 
+ K 12 (|r-r'|)Ap(r)Ap(r') + --- 



(A2) 



where Ap(r) = p(r) — p, and 



K (\t-t'\) = K(p,p;\r-r'\), 
Kl (\v-v'\) = ^lPW'P( r, );|r-r'|) 



^ii(|r-r'|) 



^I2(ir-r'|) 



dp(r) 

a 2 K(p(r),p(r');|r-r'| 



dp 2 (r) 
d 2 K( P (r), P (r>y,\r-r>\ 



(A3) 



dp{r)p{r') 

Then Eq. ([Al|) can be evaluated as a sum of separate convolutions with kernels K Q , K\, etc. 
WGC also demonstrated that only a few terms of the expansion in Eq . ( |A2l) are necessary 
to evaluate the convolution accurately for physical systems. 

The present real space method is clearly applicable to such kinetic energy functionals, 
provided K (q), Ki(q), etc., where q = k/{2k F ) and k F = (37r 2 p) 1//3 , can be fit well by 
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functions of the form of Eg. fl20|) . For the kernels other than Kq, a fitting form slightly 
different than Eq.([20|) is necessary, because for q — > oo, K (q) approaches a constant, as 
does the rational form Eq. (|20|), but the higher order kernels decay as q~ 2 . Thus they are 
fit with the modified form: 

AT 2 g2 + ... + i V 2(m _ l) ^( m -l) 



D + D 2 f + ■■■ + D 2m <?™ ■ 
Functions of this form decompose into the following terms: 



(A4) 



m „ 

J=l * J 

i.e. just as in Eq. fl2lD , but without the q 2 in the numerator. Now the Vj(r) are obtained 
by solving a modified form of Eq. (]24f) : 

1 



V 2 + S 3 



V J (r)=R J f(r), (A6) 



(2k F ) 

which is still a complex Helmholtz equation. The kinetic energy functionals of WGC have 
three parameters, a, f3, and 7. Presently, only the case of {a, (3} = |± g , 7 = 2.7, (suggested 
by favorable numerical tests) is considered. In Fig. [|, the best fit to these kernels with a 
rational function of order eight (m = 4) is shown. The quality of the fit is excellent for K Q (q) 
and Ki(q), and reasonable for the second order kernels Kn(q), K^q)- WGC have shown 
that second order terms contribute much smaller parts of Tx[p\ (which is already a small 
part of the total energy) than the zeroth and first order terms, and thus a higher fraction 
of error can be tolerated in these terms. The Pj, Qj of the decompositions of the fits to Kq 
(Eq.©) and the R jf Sj of the fits to K x , K u , and K 12 (Eq.flAlD) are given in Table glj. 
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FIG. 1: Kinetic energy functional kernel K(x), where x = multiplied by x 2 , showing 

long-ranged nature. 
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FIG. 2: The difference between the rational fitting functions K(q) (Eq.(|20|)) and K(q) (Eq. 
with the number of terms ranging from m = 1 to m = 4. The inset shows K{q). 
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FIG. 3: The fractional deviation of T^[p\ from Tk[p] for densities p(r) obtained from bulk Al 
calculations (see text). 
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FIG. 4: The fitting of the various density-dependent kernels of WGC]9|] by the forms of Eqs. (|2[ 
and ( |A4| ) with m = 4. Notice the vastly increased scale for K\, K\\, and K\2-, necessary to show 
their features, since they are negligible on the scale of Kq. 



TABLES 



TABLE I: Optimized fitting parameters for K(q) with an order-eight (m = 4) rational function 



of Eq.(21). The parameters with even indices j = 2,4 are complex conjugates of the ones given: 



P2 = P*, Pi = P3, Q2 = Q*, and Q 4 = Qt- 





i = i 


i = 3 


Pi 


0.026696 -HO. 145493 


-0.826696 + iO.691930 


Qj 


-0.818245 - i0.370856 


0.343051 - i0.689646 
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TABLE II: Lattice constant oq, in A, and bulk modulus B, in GPa, for bulk fee aluminum with the 
WGC kinetic energy functional |9| with a density-independent kernel, and with {a, (3} = | ± 
compared to those values determined with the present real-space method with m = 2, 3, and 4. 
Also shown is ag and B from the Kohn-Sham (KS) calculation. 





KS 


reciprocal space 


m = 4 


m = 3 


m = 2 


ao(A) 


4.027 


4.035 


4.030 


4.045 


4.063 


B (GPa) 


68.5 


71.9 


72.3 


68.6 


68.0 



TABLE III: Optimized fitting parameters Pj, Qj, Rj, and Sj of Eqs. (21) and (A5) for fits to the 
kernels Ko(q), K\(q), Kn(q), and K\2(q) of the WGC density-dependent kinetic energy functional 
with {a,/5} = | ± and 7 = 2.7. The parameters with even indices, j = 2,4, are complex 
conjugates of the ones given: X2 = X±, and X4 = X|, where X = P, Q, R, or <S. 





J = l 


i = 3 




Pi 


0.108403 + i0.079657 


-0.908403 + i0.439708 




Qj 


-0.470923 - i0.465392 


0.066051 - i0.259678 




Rj 


-0.030515 + i0.015027 


0.028915 - i0.008817 




Sj 


-0.597793 - i0.294130 


-0.087917 - i0.164937 




Rj 


0.008907 - i0.032841 


-0.034974 + i0.009116 


P 2 K U 


Sj 


-0.537986 - i0.233840 


-0.041565 -»0. 196662 




Rj 


0.012423 - i0.034421 


-0.031907 + i0.007392 


P 2 K 12 


Sj 


-0.511699 - i0.266195 


-0.034031 -iO.188927 
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